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ABSTRACT 

We analyze a suite of global magnetohydrodynamic (MHD) accretion disk simulations in order to 
determine whether scaling laws for turbulence driven by the magnetorotational instability, discovered 
via local shearing box studies, are globally robust. The simulations model geometrically-thin disks 
with zero net magnetic flux and no explicit resistivity or viscosity. We show that the local Maxwell 
stress is correlated with the self-generated local vertical magnetic field in a manner that is similar to 
that found in local simulations. Moreover, local patches of vertical field are strong enough to stimulate 
and control the strength of angular momentum transport across much of the disk. We demonstrate 
the importance of magnetic linkages (through the low-density corona) between different regions of the 
disk in determining the local field, and suggest a new convergence requirement for global simulations 
- the vertical extent of the corona must be fully captured and resolved. Finally, we examine the 
temporal convergence of the average stress, and show that an initial long-term secular drift in the 
local flux-stress relation dies away on a time scale that is consistent with turbulent mixing of the 
initial magnetic field. 

Subject headings: accretion, accretion disks — instabilities — MHD — turbulence 



1. INTRODUCTION 

The modern theory of accretion disks has been dom- 
inated by the discovery that angular momentum trans- 
port can be mediated by magnetohydrodynamic (MHD) 
turbul ence driven by the magnetor otational instability 
(MRI: iBalbus k Hawlevlll99lL Il998l ). Although analytic 
treatments of the MRI suffice to establish that the in- 
stability exists and rapidly develops toward turbulence, 
numerical work has been at the forefront of the effort to 
characterize the resulting angular momentum transport. 
Simulations of the MRI require difficult compromises be- 
cause arguably important scales span a wide range be- 
tween the global scale A ~ r, an intermediate scale A ~ h 
(where h is the disk scale height) that roughly defines the 
scale over which shear dominates turbulent fluctuations, 
and viscous and resistive dissipation scales A„ and A,,. 
As in many other astrophysical problems, the latter are 
usually so small that in the disks under consideration it 
is impossible to run simulations that capture the physical 
dissipative scales. 

To date, the bulk of our numerical understanding of 
the MRI has been derived from local shearing-box sim - 
ulations (|Hawlev et all 119951 : [Brandenburg et ; all Il995f ). 
In a local shearing-box model, one studies the local dy- 
namics of an orbiting patch of the accretion flow using 
a (co-rotating) Cartesian coordinate system by including 
Coriolis forces and shearing boundary conditions. Since 
only a small patch of the disk is being simulated, this 
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approach maximizes the separation between the interme- 
diate driving scale of the turbulence and the dissipative 
scales. For our purposes, the result that in local sim- 
ulations the saturation level of magnetic fields and the 
strength of angular momentum transport are found to 
scale with the vertical flux threading the simulation do- 
main (jHawlev et alJl!995t iSano et al.N 2004: Pc ssah et all 
I2007t ) will be of particular interest. 

Local simulations h ave know n lim itations 
(|Reeev k Umurhanl [20081 : iBodo et all [2001 . By 
construction, they enforce the local conservation of 
quantities (in particular the net magnetic flux) that 
in reality are only globally conserved. In many older 
implementations they also enforce periodicity (in 
radius, azimuth, and in some instances also height) 
on a scale that may be small enough to impact the 
result s. Recently, a number of auth ors dDavis et al.l 
120091: iGuan et alj|2009t Uohansen et al.ll2009D have used 
larger than usual shearing-box simulations to quantify 
whether these limitations matter for practical purposes. 
In this paper we address the same problem from the 
other direction. We analyze small patches of global disk 
simulations in an attempt to determine whether the 
disk behaves as if it were a collection of shearing-boxes. 
Our specific goal is to ascertain whether the relationship 
between local r</>-component of the magnetic stress and 
vertical magnetic flu x that is found in local simulations 
(jPessah et alj 120071 ) is recovered in local patches of 
global simulations. As a result, we uncover the impor- 
tance of the magnetic connectivity of the disk and the 
need to fully capture the vertical extent of the corona. 

This paper is organized as follows. §2 briefly describes 
the global simulations that we employ. §3.1 describes 
the vertical structure of the disks we simulate with an 
emphasis on the distribution of magnetic flux. We dis- 
cuss, in Section §3.2, the results of our study of instanta- 
neous correlations between magnetic flux and stress. The 
flux-stress connection is explored in more detail in Sec- 



2 



tion §3.3 (where we discuss the nature of the transition 
point in the flux-stress relation) and Section §3.4 (where 
we examine the temporal correlation of stress and ver- 
tical flux in co-moving patches). Sections §3.5 and §3.6 
describe the dependence of angular momentum transport 
on vertical domain size and time respectively. Section §4 
presents our conclusions. 

2. GLOBAL SIMULATIONS OF THIN ACCRETION DISK 

The sim ulations used here are very similar to those de- 
scribed bv [Reynolds fc Miller! (l2009t). We use the ZE US- 
MP code (|Stone fc Normanl 119921 : lHaves et all 120061 ) to 
solve the equations of ideal MHD in three dimensions. 
We modified the basic version of the code to incorporate 
a Paczynski-Wiita pseudo-Newtonian gravitational po- 
tential (as a first approximation of th e gravitational field 
abou t a Schwarzschild black hole; iPaczynski fc Wiital 
1980) and performed the simulation in cylindrical po- 
lar coordinates. Our simulations are ideal MHD in the 
sense that no explicit resistive or viscous dissipation is in- 
cluded; all dissipative processes are due to the discretiza- 
tion of the spatial domain and hence occur close to the 
grid scale. Furthermore, we integrate an internal energy 
equation assuming an adiabatic equation of state with 
7 = 5/3. Energy is lost from the domain when magnetic 
fields undergo numerical reconnection. 

The initial disk is in a state of Keplerian rotation (with 
respect to the pseudo-Newtonian potential), and is in 
vertical hydrostatic equilibrium with a constant scale- 
height h. Thus, the initial density, pressure and velocity 
field is 
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where r represents the cylindrical radius, R = \Jr 2 + z 2 , 
r g = GM/c 2 and h = 0.05r isco = 0.3r 9 . We set the 
initial midplane density to be po(r) = 1 beyond the in- 
nermost stable circular orbit (ISCO) at r isco = 6r s , and 
Po{r) — within the ISCO. The radial boundary condi- 
tions correspond to zero-gradient outflow, while periodic 
boundaries were imposed on both the <f> and z-boundaries 
(this last choice, which is made in order to avoid field- 
line snap ping and other numerical problems, is discussed 
further in iReynolds & Milled I2009T ). 

The initial magnetic field is specified in terms of a vec- 
tor potential of the form, 
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Here / is an envelope function that is unity in the disk 
body and smoothly goes to zero away from the main body 
so as to avoid unphysical interactions with the bound- 
aries. This results in a magnetic field topology consist- 
ing of distinct poloidal field loops of alternating orien- 
tation throughout the main body of the simulation. Of 
importance to the current discussion is that there is no 



net vertical magnetic flux threading the disk as a whole. 
The constant Aq is chosen to ensure that the magnetic 
field strength is normalized so that the ratio of volume- 
integrated gas and magnetic pressure (3 10 3 . 

A set of simulations were run to span a range of nu- 
merical parameters, specifically varying the vertical and 
radial resolution, as well as the vertical and azimuthal 
extent of the domain. A comparison of simulations with 
varying azimuthal extents suggest a negligible depen- 
dence on this parameter, and as such all the simula- 
tions considered here use the same 30° wedge-shaped az- 
imuthal domain. The vertical domain size was found to 
be important, and we will discuss the role of this pa- 
rameter later in this paper. A detailed study of the de- 
pendence on resolution, requiring much greater compu- 
tational expense, is deferred to a later work. All simula- 
tions presented here have a radial domain r £ (4r fl , 16r g ). 
Details of the simulations considered are given below in 
Table 1. 

3. RESULTS 
3.1. Vertical Structure of Thin Disks 

Our primary goal is to study the instantaneous corre- 
lation between stress and magnetic flux within the simu- 
lated disk. In a stratified disk this correlation may vary 
with height above the disk midplane, and hence we start 
by considering how the mean magnetic field structure 
varies vertically in our simulations. 

To analyze the simulations the principle quantity of 
interest is the r — <fi component of the Maxwell stress 
tensor, 
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whic h dominates MRI-driven angular momentum trans- 
port (Bal bus fc Hawlevlll998l ). In a turbulent disk, both 
M r( j, and other physical quantities of interest are com- 
plicated functions of space and time. To make sense 
of them we use temporal and spatial averages. Run 
Thin.M-Res . 6z has the longest duration of any of our 
simulations, and we use this run to construct a rep- 
resentative vertical profile of the magnetic structure of 
the disk. To reduce the effects of spatial intermittency 
in the turbulence, we azimuthally average over the en- 
tire domain and average over a small radial range cen- 
tered about a fiducial radius in the body of the disk ( 
8r g — h < r < 8r g + h) . To smooth out the temporal 
variability and isolate the behavior of the disk in a sat- 
urated turbulent state we time average over 400 ISCO 
orbits starting at orbit 50. 

The results are given in Figure [1] which shows the ver- 
tical profiles of the relevant quantities scaled to their 
maximum values. Our interest in vertical structure is 
predominantly in the magnetic fields, and in particular 
the vertical magnetic flux, B z , and the magnetic stress, 
M rc f,. However, we also plot the density, p, to highlight 
the contrast between the relatively unmagnetized mid- 
plane of the disk and the sparse magnetized "corona" 
away from the midplane. The obvious reason for the for- 
mation of these two disparate regions is magnetic buoy- 
ancy resulting from the effect of vertical gravity, but this 
may be overly simplistic. Also of interest is the double- 
peak vertical profile of vertical flux and stress. Broadly 
similar results are seen in a subset of the stratified lo- 
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Run ID 


Resolution 


Vertical Extent 


Total Orbits 


Mean olm 


Standard Deviation 




(R,0,z) 


( in r g ) 


(at T*isco) 




of a M 


Thin.M-Rcs.12z 


(240,32,512) 


12 


122 


0.008 


0.0011 


Thin.M-Res.6z 


(240,32,256) 


6 


664 


0.0086 


0.0016 


Thin.M-Res.3z 


(240,32,128) 


3 


112 


0.0061 


0.001 



TABLE 1 

The spatial resolution, vertical domain size, duration of the simulations analyzed in this work. Also included are the 
mean and standard deviation of «m, defined by equation 1141 between 50 and 100 orbits at r isco . 
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Fig. 1. — Vertical structure of physical quantities averaged over 
four hundred orbits for run Thin.M-Res . 6z. 




kt 4 i<r 3 icr 2 io"' io° 



FlG. 2. — Flux versus stress averaged over four hundred orbits 
for run Thin.M-Res.6z. 



cal simulations of lMiller fc Stone] (12000D and a subs et of 
the global simulations of iFromang fc Nelson! (|2006l ). It 
is interesting that the vertical location of the peak field 
and stress seem to approach constant values rather than 
growing monotonically with time. Whether the region 
of strong flux is trapped, possibly due to magnetic ten- 
sion from field lines connecting it to the midplane, or 
is continually dissipating while outflowing, is currently 
unclear. 

The local flux-stress relationship we will compare our 
results to is based on unstratified shearing box simu- 
lations and the rich vertical magnetic structure due to 
stratification in our global disk means that there is nec- 
essarily some ambiguity in the comparison. In what 
follows, we therefore analyze the flux-stress relationship 
within the global simulation not just at the midplane but 
also as a function of height. One should note that the off- 
set of the peak vertical flux and stress from z = means 
that high values of flux and stress are only accessible 
away from the midplane. 

3.2. A local flux-stress correlation in a global disk 

Although our simulations have zero net vertical flux, 
small patches of the disk are instantaneously threaded 
by a vertical field. We seek to determine whether the 
Maxwell stress tracks this transient vertical field in the 
same way as it would in a local simulation where the ver- 
tical field is persistent (jHawlev et alj|1995t IPessah et al.l 
I2007T ) . To proceed, we break up the global simulation do- 
main at each timestep into several hundred small cylin- 
drical wedges of size Az = Ar = h and A0 « 0.1. Within 
each wedge we average to obtain a single estimate of 
the magnetic stress (M r( /,/p, normalized to the local gas 
pressure) and the local vertical flux, which we express in 



terms of the wavelength of the most unstable MM mode, 

where va z = Bz/^/Airp is the Alfven speed correspond- 
ing to the vertical field component. Note that because 
we scale the stress by the local pressure (a decreasing 
function of height), we immediately introduce a height 
dependence. Choosing instead to scale by the midplane 
pressure (only a function of radius) still yields a height 
dependence and as such we are confident that the height 
dependence seen is not solely a consequence of the pres- 
sure scaling. 

To avoid early transients, our analysis excludes the pe- 
riod prior to the first 50 ISCO orbital periods. To im- 
prove our statistics (and to give a measure of the conver- 
gence of our results given the finite duration of the run) 
we consider wedges that are centered at z = ±[0, 1, 2, 3]h 
and plot results for the samples separately. The result- 
ing pairs of flux-stress values from all of the wedges and 
all of the snapshots in time were binned according to 
(logarithmic) vertical flux in order to diagnose trends. 

The resulting flux-stress relations for all the simula- 
tions considered are broadly similar. Our best statistics 
come from the long duration run Thin.M-Res . 6z, and 
the flux-stress relation for this case is plotted in Figure [2j 
The stress is observed to be flat for weak vertical fields 
(small Amri), while for larger field strengths we have ap- 
proximately M r( p oc Amri- This may be compared with 
the local scaling relation derived from unst ratified simu- 
lations, which can be written in the form (IPessah et al.l 
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where A and L are the vertical grid cell size and total box 
size respectively. The pressure scale height H is given by 
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where Slo is the Keplerian angular velocity and c s the 
sound speed. Note the distinction between the locally 
defined H and the globally constant h. In comparing 
our results to those obtained for local simulations, we 
consider L to be the size of the wedge and H to be the 
locally defined pressure scale height. 

Our simulations do not spontaneously develop vertical 
fields strong enough to quench the MRI (and hence we 
do not sample the Amri > L regime), but the behavior 
of patches at low and intermediate vertical fluxes is qual- 
itatively the same as that found in local si mulations. In 
contra st to the local (unstratified) results of lPessah et al.l 
(|2007f ) is the strong dependence on height in our strati- 
fied simulations. In addition to the fact that the largest 
values of flux are only accessible at large heights is the 
fact that the stress response to flux is also height de- 
pendent. There is a height-dependence of the transition 
point between low and intermediate flux, as well as in the 
slope of the intermediate flux regime. Of particular note 
is the location of the transition point itself. The vertical 
line marked in Figure [2] is given by Amri = A/20, and 
approximately marks the location of the transition point 
at z = ±H. This stands in contrast to the transition 
point for local unstratified simulations, Amri = A z . We 
return to a discussion of the physics of this transition in 
Section 13.31 

The fact that transient self-generated vertical flux is 
able to stimulate the local stress in the same manner as 
occurs in local models is primarily a formal result, al- 
though it does lend some credence to models in which 
patches of vertical field are assum ed to have a physical 
identity ([Spruit fc Uzdenskvll2005() . Of greater import is 
the observation that, across much of the disk, the self- 
generated field is strong enough to fall into the linear 
regime of the flux-stress relation. Figure [3] shows the 
distribution of flux through patches in all three. The 
vertical flux distributions are approximately symmetric 
in log-space. We find that, at any instant in time, about 
half of the area of the disk is threaded by a field strong 
enough to control the stress. We interpret this to mean 
that in a zero-net field global simulation, much of the 
disk sees a field strong enough to control its dynamics. 
In other words, the dynamics of the disk is strongly influ- 
enced by the connectivity of the self-generated magnetic 
field between different patches of the disk. 

In these results we again see the importance of the 
magnetized region away from the midplane. Not only is 
this the location of the largest values of flux, and thus 
stress, but for the same flux the stress response is higher. 
This increased stress response to flux suggests that the 
corona is important not just as a warehouse of magnetic 
energy, but also has the ability to use this magnetic en- 




FlG. 3. — Normalized flux distributions for all simulations, com- 
puted between 50-100 ISCO orbits. The dashed vertical line in- 
dicates, approximately, the strength of vertical field above which 
magnetic stress in the disk is stimulated. 

ergy more efficiently to induce angular momentum trans- 
port. One possible explanation for this efficiency is the 
ability of the corona to mediate magnetic links through 
radially disparate regions of the disk due to the presence 
of coronal field loops. In Section [3~5l we consider the ef- 
fect of vertical domain size on the saturated turbulent 
stress as a means of diagnosing whether truncating the 
corona affects the dynamics of a disk. Before that, how- 
ever, we discuss the nature of the transition point in the 
flux-stress relation as well as study the connection be- 
tween the local vertical magnetic fields and stresses via 
their temporal correlation. 

3.3. The nature of the transition point in the flux-stress 

relation 

To reiterate, the flux-stress relation obtained from our 
global simulation shows a transition at approximately 
Xmri ~ A/20, in contrast with Xmrt ~ A foun d 
from local unstratified simulations (Pes sah et alJ 12007). 
Neither of these transitions can correspond straightfor- 
wardly to the condition that the fastest growing MRI 
mode is resolved. To resolve a mode in a ZEUS-likc 
scheme requires that the wavelength is spanned by at 
least ~ 8 computational zones. Since Xmri is, by con- 
struction, the wavelength of the fastest growing mode 
corresponding to the net vertical magnetic field, the con- 
dition that the fastest growing MRI mode is resolved cor- 
responds to Xmri ~ 8A. Our result then implies that 
magnetic fields that are very weak - in the sense that 
their fastest growing mode cannot be resolved - nonethe- 
less have an important influence on the dynamics of our 
global disk. 

We do not have a quantitative explanation of why 
^mri ~ A/20. On general grounds, however, we note 
that we would expect that the transition point would lie 
at Xmri *C 8A. A given vertical field is unstable not just 
to the fastest growing mode, but also to a whole spectrum 
of slower-growing modes that have longer wavelengths 
that are more easily resolvable numerically. Plausibly, 
the transition point will then correspond to the condi- 
tion that we resolve the slowest growing mode that grows 
appreciably before it is truncated by non-linear coupling 
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to other MRI modes or some other aspect of the physics 
(e.g. a dynamo cycle). If this is the case, then it is un- 
surprising that the transition point varies between local 
and global simulations, since the time scale available for 
a mode to grow may well depend on the presence or ab- 
sence of a low density disk corona within which the MRI 
is not active. 

To consider this more quantitatively, consider purely 
vertical MRI modes (k r — h<p — 0) in a thin accretion 
disk (so that radial gradients of pressure and entropy 
can be neglected). Let the spacetime dependence of the 
modes be e 'M- fcz ). The disp ersion relation for these 
modes (Balb us fc Hawlevlll99lD reads 
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where Q 2 = uj 2 — k 2 v 2 A and k is the radial epicyclic (an- 
gular) frequency. We wish to examine modes with wave- 
lengths much longer than the fastest growing mode, i.e., 
with Ifcv^/fiol 1- Rewriting in terms of the growth 
rate, a = —iuj and expanding the dispersion relation to 
lowest order in k 2 v 2 A /Vi^ gives 
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Suppose that a given mode can grow exponentially for 
a time r before it is truncated by mode coupling or 
some other unspecified physical process. Then, the slow- 
est growing mode that actually experiences significant 
growth (hereafter, the slowest appreciably growing mode 
[SAGM]) has (T sagrn — 2tt/t and a wavenumber given by 
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Our hypothesis is that the transition point in the flux- 
stress relation corresponds to the point where the slowest 
appreciably growing mode is just resolvable, i.e., where 
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where r' = r/t or f,, t or i, being the orbital period at that 
radius. 

Equation [T3] offers some insight into the transition 
point found in the flux-stress relations that we have 
be en considerin g ,. In the local unstratified simulations 
of iPessah et al.l^2007t) . the implicit potential is Newto- 
nian (k 2 = Qq) and we find that Xmri ~ A, imply- 
ing r ~ t or b. In the global simulations presented here, 
we find the transition point at Xmri ~ A/20 which 
(accounting for the fact that k 2 < fig in the pseudo- 
Newtonian potential) gives r ~ 5 — 10t or t,. Thus, within 
the framework of this argument, the difference in the lo- 
cation of the transition point between the local unstrat- 
ified and the global simulations is due to a difference in 
the robustness of the long wavelength and slowly grow- 
ing modes; slowing growing modes appear to be able to 
grow for longer within the global simulation before being 
truncated. The nature of this difference, which must be 
closely related to the saturation of the turbulent state, 



is beyond the scope of this paper and will be explored in 
future work. 

3.4. Temporal correlations in flux-stress 

The results of SJ3T2] suggests that the (fluctuating) mag- 
netic flux threading a local patch of the disk determines 
the r — 4> component of the magnetic stress generated 
by the turbulence in that patch. If the vertical mag- 
netic flux is indeed the causal agent in determining the 
stress, we expect a temporal lag between fluctuations 
in the magnetic flux and the resulting variations in the 
stress. On the basis of experiments with local simulations 
(jHawley et al.l fl996h . we expect this lag to be approxi- 
mately two (local) orbital periods. Thus, we expect the 
temporal lag to increase with radius in the disk due to 
the increasing orbital period. 

To search for this lag, we use Thin.M-Res . 6z and out- 
put the 3-d structure of the disk once every 0.1 ISCO or- 
bits during the interval between 50-90 ISCO orbits (this 
is 10 times the nominal data output rate). Using these 
400 snapshots of the disk structure, we then computed 
the instantaneous vertical magnetic fluxes and magnetic 
r — <p stresses in families of co-moving wedges at three 
radii r £ {8r g , 10r g , 12r g }. The azimuthally averaged 
value of i>0 at each radius was used to track a given co- 
moving wedge between timesteps. This procedure is not 
fully lagrangian, because it does not account for the ra- 
dial movement or fluctuating azimuthal velocity of a co- 
moving patch, but we expect these effects to be negligible 
for the short timeframe under consideration. The time- 
series of magnetic flux and stress for each wedge were 
then cross-correlated and, finally, the cross-correlations 
for all wedges at a given radius were averaged. 

The resulting averaged temporal cross-correlations are 
shown in Fig. [4] At each radius we see a strong instan- 
taneous correlation, likely due to the immediate shear- 
ing of perturbed vertical fields. However, in general, the 
cross-correlation is biased toward positive lag. This is 
consistent with what we would expect, namely that the 
presence of vertical flux will feed the MRI and result in 
enhanced transport. Of note is the fact that the inner- 
most radius considered, R = 8r g exhibits a double peak 
structure whereas this is unresolved at higher radii. Also 
peculiar is the fact that the outer-most radius, R = 12r g 
is significantly less biased towards positive lag than the 
other radii under consideration. A further exploration of 
these issues is beyond the scope of this paper, and will 
be explored in future work employing orbital advection 
algorithms and test-particle tracers in order to correctly 
follow the evolution of a local patch. 

3.5. The dependence of stress on vertical domain size 

When considering the correlation between vertical 
flux and stress we work directly with M r( p, but when 
studying possible trends in average stress with do- 
main size we instead defi ne an effective a-parameter 
( Shak ura fc Sunvaevj|1973l) . For each snapshot we define 
a density-weighted spatial average via, 
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Fig. 4. — Cross-correlation of flux and stress evaluated in co- 
moving patches of the disk. Positive lags imply that changes in 
the vertical flux occur, on average, prior to changes in the magnetic 
stress. 

The restriction on the radial range of the averaging is 
designed to ignore the plunging region of the accretion 
flow (r < 6r g ) and any effects related to the outer radial 
boundary. Density weighting is used in the vertical di- 
rection to take into account the low density, highly mag- 
netized regions while still allowing the dominant contri- 
bution to the integral to come from the denser mid-plane 
of the disk. 

A comparison of otyi and its dependence on vertical 
domain is given in Figure [5] The initial growth phase 
of the MRI is unaffected by the vertical domain as ex- 
pected, since all the simulations considered have the same 
vertical resolution and can thus resolve the same unsta- 
ble MRI modes. Over the course of the simulations, the 
larger vertical extent simulations have, in general, larger 
values of om- We attribute this effect to stifling of the 
growth of the magnetized regions in the smaller verti- 
cal extent simulation. However, the long-term effects of 
vertical extent are ambiguous. Whether the simulations 
converge to the same otyi or the apparent convergence 
is a result of short-term variability is unclear from the 
current simulations. Longer simulations will need to be 
carried out to determine the vertical domain size that 
is needed in order to reliably capture the dynamics of a 
global disk. 

3.6. Long term behavior 

In the one case of run Thin.M-Res . 6z, the disk was 
simulated for 664 ISCO orbits. This simulation allows 
us to search for long-term trends in the dynamics of the 
disk. As shown in Figure [5] there is a slight downward 
drift in olm over time. The same temporal trend is also 
evident in the flux-stress relationship. Figure [7] shows 
the flux-stress relationship averaged in 100 orbit blocks 
starting at 50 ISCO orbits. During the first 300 orbits, 
there appears to be a secular drift in the flux-stress rela- 
tion. The linear (high-flux) part of the relation achieves 
a steady state relatively quickly (only the first time block 
between 50 and 150 ISCO orbits shows significant differ- 
ences), but the flat (low-flux) part of the relation con- 
tinues to fall until it too achieves a steady state at ap- 
proximately 350 ISCO orbits into the run. Associated 
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Fig. 5. — Behavior of om an< i its dependence on vertical domain. 

with this, the "knee" in the flux-stress relation appears 
to move to smaller fluxes. 

In essence, this result says that low-flux regions still 
support (small) stresses at early times but that those 
stresses decay over a period of several hundred ISCO or- 
bits. We ascribe this to stresses associated with a sheared 
residual of the initial magnetic field configuration which 
are "mixed away" on a relatively long timescale. Our 
initial field configuration threads the midplane with re- 
gions of net magnetic flux which alternate with a radial 
periodicity of 5h. Radial Fourier transforms of the mid- 
plane azimuthally-averaged B z do indeed find a (weak) 
periodicity corresponding to the initial field even once 
the turbulence is fully developed. This periodic compo- 
nent grows weaker and is no longer detectable at approx- 
imately the same time that the flux-stress curve achieves 
steady-state. These observations further suggest that 
residual flux from the initial conditions is responsible for 
the long term variability. 

Assuming that a long-lived residual of the initial mag- 
netic field is the driving mechanism for this phenomenon, 
we can recover the time required to achieve the steady 
state from elementary arguments. The time needed to 
turbulently diffuse together two patches of oppositely di- 
rected flux separated by a radial distance Ar = 2.5h is 
given by 

Ar 2 

^mix ~ i (15) 

VeS 

where rj e s is the effective turbulent resistivity. If we de- 
fine Pr miC ff as the effective turbulent magnetic Prandtl 
number (i.e. the ratio of the effective turbulent viscosity 
to the effective turbulent resistivity), we can write 

?7eff = Pr~* cS a M c s h, (16) 

where c s is the sound speed. We can then write the 
mixing time as 

P7Veff/Ar\ 2 

* mbc 2^{-h) toTh > (17) 

where t or b is the local orbital period and we have 
used the fa ct that h/c s ~ T.l v * . ~ ^orb/27r. Using 
fVm.eff = 1 (|Guan fc Gammiel 120091 : iLesur fc Longarettil 
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Fig. 6. — Long term behavior of olm for run Thin.M-Res . 6z. 




Fig. 7. — Flux versus stress averaged over 100 orbit blocks for 
run Thin.M-Res . 6z. Vertical centering of z = h. 



[20091 iFromang fc Stondl2009l ) and a M = 0.005 suggests 
that the memory of the initial conditions will be lost on 
a timescale of t m i x ~ 200£ O rt>- This crude estimate is in 
reasonable agreement with the timescale on which we see 
the flux-stress relationship achieve a stationary state. 

4. CONCLUSIONS 

It has been a long-held ansatz that one can extract 
and model the dynamics of a local patch of an accre- 
tion disk and obtain results (for the angular momentum 
transport, for example) that have meaning for the disk 
as a whole. By examining local patches of a high reso- 
lution global disk simulation, we have provided a direct 
test of this notion. We have shown that MRI-driven tur- 
bulence in global geometrically thin accretion disks be- 
haves in a way consistent with scaling laws derived for 
local simulations. In particular, we find that global disks 
display a local flux-stress relation qualitatively similar 
to that found in l ocal simulations (|Hawlev et al.l H9951 ; 
iPessah e~ al. 2007). However, other aspects of the global 
models are distinctly different: 



1. Even though we model a global accretion disk that 
has zero net magnetic field, any given patch of the 
disk is threaded by a net magnetic flux resulting 
from the self-generated field in the MRI-dynamo. 
Across much of the disk, the local flux is strong 
enough to have a controlling effect on the local 
stress. Thus, our zero-net field global disk is be- 
having as a collection of net field local patches. 

2. The normalization, slope, and location of the 
"knee" of the flux-stress relationship changes with 
vertical height in the accretion disk. This amplifies 
the role of the off-midplane region (h < \z\ < 2h) 
of the disk; not only does this region have stronger 
vertical magnetic fields than the midplane, but 
a given vertical field induces stronger magnetic 
stresses. The result is a strong enhancement of 
magnetic stress well off the midplane of the disk. 

3. The transition point (or "knee") in the flux-stress 
relation occurs as significantly smaller fluxes in the 
global simulation as compared to the local unstrat- 
ified simulations. We relate this transition point to 
the ability of the simulation to marginally resolved 
the slowest appreciably growing mode. 

4. Angular momentum transport (i.e. olm) in the 
global disk appears to be impeded if the vertical do- 
main size of the simulation is too small; we found 
significant differences between our z = ±5/i and 
z = ±10h cases. On the other hand, the z — ±10h 
and z — ±20/i cases appear very similar suggesting 
convergence has been achieved. Given that mag- 
netic linkages between different patches of the disk 
appear to be crucial for determining the local flux 
(and hence the local stress), and that such link- 
ages are made through the low-density corona of 
the disk, such a sensitivity to the vertical domain 
size is not surprising. 

5. Analysis of our long simulation (which ran for 664 
ISCO orbits) reveals long-term secular trends. In 
particular, there is a secular drift of the flux-stress 
relationship such that the knee of the relationship 
moves to smaller fluxes and the low-flux normal- 
ization decreases. This secular drift stabilizes after 
approximately 300 ISCO orbits. We attribute this 
to stresses associated with a long-lived residual of 
the initial magnetic field configuration. 
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